Higher-order numerical methods for stochastic simulation of chemical reaction systems

In this paper, using the framework of extrapolation, we present an approach for obtaining higher-order τ-leap methods for the Monte Carlo simulation of stochastic chemical kinetics. Specifically, Richardson extrapolation is applied to the expectations of functionals obtained by a fixed-step τ-leap algorithm. We prove that this procedure gives rise to second-order approximations for the first two moments obtained by the chemical master equation for zeroth-and first-order chemical systems. Numerical simulations verify that this is also the case for higher-order chemical systems of biological importance. This approach, as in the case of ordinary and stochastic differential equations, can be repeated to obtain even higher-order approximations. We illustrate the results of a second extrapolation on two systems. The biggest barrier for observing higher-order convergence is the Monte Carlo error; we discuss different strategies for reducing it.